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1 Incorporation By Reference 

2 U. S. Patent Application Serial No. 08/880,418, filed June 23, 1997, in the names of Martin 

3 Grabenstein, et al., entitled "System And Method For Generating Pixel Values For Pixels In An 

4 Image Using Strictly Deterministic Methodologies For Generating Sample Points," (hereinafter 

5 referred to as the Grabenstein application) assigned to the assignee of this application, incorporated 

6 by reference. 

**% Field Of The Invention 

fn 
o5 

j| The invention relates generally to the field of computer graphics, and more particularly to 

fft 

21 systems and methods for generating pixel values for pixels in the image. 

y a 

r™ 

tQ Background Of The Invention 

P| In computer graphics, a computer is used to generate digital data that represents the 

P projection of surfaces of objects in, for example, a three-dimensional scene, illuminated by one or 

1 3 more light sources, onto a two-dimensional image plane, to simulate the recording of the scene by, 

14 for example, a camera. The camera may include a lens for projecting the image of the scene onto 

15 the image plane, or it may comprise a pinhole camera in which case no lens is used. The two- 

1 6 dimensional image is in the form of an array of picture elements (which are variable termed "pixels" 

17 or "pels"), and the digital data generated for each pixel represents the color and luminance of the 

1 8 scene as projected onto the image plane at the point of the respective pixel in the image plane. The 

19 surfaces of the objects may have any of a number of characteristics, including shape, color, 

20 specularity, texture, and so forth, which are preferably rendered in the image as closely as possible, 

21 to provide a realistic-looking image. 

22 Generally, the contributions of the light reflected from the various points in the scene to the 

23 pixel value representing the color and intensity of a particular pixel are expressed in the form of the 

24 one or more integrals of relatively complicated functions. Since the integrals used in computer 
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graphics generally will not have a closed-form solution, numerical methods must be used to evaluate 
them and thereby generate the pixel value. Typically, a conventional "Monte Carlo" method has been 
used in computer graphics to numerically evaluate the integrals. Generally, in the Monte Carlo 
method, to evaluate an integral 

1 

(/) = \f(x)dx (1) 

0 

where f(x) is a real function on the real numerical interval from zero to 1, inclusive, first a number 
"N" statistically-independent random numbers xi , i=l, N, are generated over the interval. The 
random numbers xi are used as sample points for which sample values f(xi) are generated for the 

function f(x), and an estimate / for the integral is generated as 

</>»/ = £l; :,/(*,) » 



As the number of random numbers used in generating the sample points f(xi) increases, the value 
of the estimate f will converge toward the actual value of the integral i^f^j . Generally, the 

distribution of estimate values that will be generated for various values of "N," that is, for various 
numbers of samples, of being normal distributed around the actual value with a standard deviation 
a which can be estimated by 

if the values Xj used to generate the sample values f(x { ) are statistically independent, that is, if the 
values ^ are truly generated at random. 

Generally, it has been believed that random methodologies like the Monte Carlo method are 
necessary to ensure that undesirable artifacts, such as Moire patterns and aliasing and the like, which 
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1 are not in the scene, will not be generated in the generated image. However, several problems arise 

2 from use of the Monte Carlo method in computer graphics. First, since the sample points x i used in 

3 the Monte Carlo method are randomly distributed, they may clump in various regions over the 

4 interval over which the integral is to be evaluated. Accordingly, depending on the set of random 

5 numbers which are generated, in the Monte Carlo method for significant portions of the interval 

6 there may be no sample points x { for which sample values f(Xj) are generated. In that case, the error 

7 can become quite large. In the context of generating a pixel value in computer graphics, the pixel 

8 value that is actually generated using the Monte Carlo method may not reflect some elements which 
J£ might otherwise be reflected if the sample points x { were guaranteed to be more evenly distributed 
tti over the interval. This problem can be alleviated somewhat by dividing the interval into a plurality 
||j of sub-intervals, but it is generally difficult to determine a priori the number of sub-intervals into 

# which the interval should be divided, and, in addition, in a multi-dimensional integration region 

8P 

11 (which would actually be used in computer graphics, instead of the one-dimensional interval 

T4 described here) the partitioning of the region into sub-regions of equal size can be quite complicated. 



In addition, since the method makes use of random numbers, the error 



/ - (/)| (where |x| 



jj| represents the absolute value of the value "x") between the estimate value f and actual value 

17 is probabilistic, and, since the error values for various large values of "N" are close to normal 

18 distribution around the actual value {f ) , only sixty-eight percent of the estimate values f that 

19 might be generated are guaranteed to lie within one standard deviation of the actual value (/) . 

20 Furthermore, as is clear from equation (3), the standard deviation a decreases with increasing 

2 1 numbers of samples N, proportional to the reciprocal of square root of "N" (that is, ^/^y^* )■ Thus, 

22 if it is desired to reduce the statistical error by a factor of two, it will be necessary to increase the 

23 number of samples N by a factor of four, which, in turn, increases the computational load that is 

24 required to generate the pixel values, for each of the numerous pixels in the image. 
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1 Additionally, since the Monte Carlo method requires random numbers, an efficient 

2 mechanism for generating random numbers is needed. Generally, digital computers are provided 

3 with so-called "random number" generators, which are computer programs which can be processed 

4 to generate a set of numbers that are approximately random. Since the random number generators 

5 use deterministic techniques, the numbers that are generated are not truly random. However, the 

6 property that subsequent random numbers from a random number generator are statistically 

7 independent should be maintained by deterministic implementations of pseudo-random numbers on 

8 a computer. 

9 The Grabenstein application describes a computer graphics system and method for generating 
J§ pixel values for pixels in an image using a strictly deterministic methodology for generating sample 

Pi points, which avoids the above-described problems with the Monte Carlo method. The strictly 

on 

jj2 deterministic methodology described in the Grabenstein application provides a low-discrepancy 

p sample point sequence which ensures, a priori, that the sample points are generally more evenly 

N- distributed throughout the region over which the respective integrals are being evaluated. In one 

|§ embodiment, the sample points that are used are based on a so-called Halton sequence. See, for 

f§ example, J. H. Halton, Numerische Mathematik, Vol. 2, pp. 84-90 (1960) and W. H. Press, et al., 
Numerical Recipes in Fortran (2d Edition) page 300 (Cambridge University Press, 1992). In a 

p 

£8 Halton sequence generated for number base "p," where base "p" is a selected prime number, the "k- 

19 th" value of the sequence, represented by Hp k is generated by use of a "radical inverse" operation, 

20 that is, by 

21 (1) writing the value "k" as a numerical representation of the value in the selected base 

22 "p," thereby to provide a representation for the value as D M D M _ X • • • D 2 D X , where 

23 D m (m=l, 2, M) are the digits of the representation, 

24 (2) putting a radix point (corresponding to a decimal point for numbers written in base 

25 ten) at the least significant end of the representation D M D M _ X • • • D 2 D l written in 

26 step (1) above, and 
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1 (3) reflecting the digits around the radix point to provide 0. D x D 2 • • • D M _ x D M , which 

2 corresponds to Hp k . 

3 It will be appreciated that, regardless of the base "p" selected for the representation, for any series 

4 of values, one, two, ... "k," written in base "p," the least significant digits of the representation will 

5 change at a faster rate than the most significant digits. As a result, in the Halton sequence Hp 1 , Hp 2 , 

6 . . .Hp k , the most significant digits will change at the faster rate, so that the early values in the sequence 

7 will be generally widely distributed over the interval from zero to one, and later values in the 

8 sequence will fill in interstices among the earlier values in the sequence. Unlike the random or 
J) pseudo-random numbers used in the Monte Carlo method as described above, the values of the 

# Halton sequence are not statistically independent; on the contrary, the values of the Halton sequence 

03 

P are strictly deterministic, "maximally avoiding' 1 each other over the interval, and so they will not 

|p clump, whereas the random or pseudo-random numbers used in the Monte Carlo method may clump. 

m 

L3 It will be appreciated that the Halton sequence as described above provides a sequence of 

J4 values over the interval from zero to one, inclusive along a single dimension. A multi-dimensional 

B5 Halton sequence can be generated in a similar manner, but using a different base for each dimension. 

M A generalized Halton sequence, of which the Halton sequence described above is a special 

H case, is generated as follows. For each starting point along the numerical interval from zero to one, 

1 8 inclusive, a different Halton sequence is generated. Defining the pseudo-sum xe p y for any x and y 

1 9 over the interval from zero to one, inclusive, for any integer "p" having a value greater than two, the 

20 pseudo-sum is formed by adding the digits representing "x" and "y" in reverse order, from the most- 

2 1 significant digit to the least-significant digit, and for each addition also adding in the carry generated 

22 from the sum of next more significant digits. Thus, if "x" in base "p" is represented by 

23 O.X x X 2 ~- X M _ X X M , where each "X m M is a digit in base "p," and if "y" in base "p" is represented 

24 by 0. Y X Y 2 • • • Y N _ X Y N , where each "Y n " is a digit in base "p" (and where "M," the number of 

25 digits in the representation of "x" in base "p", and "N," the number of digits in the representation of 

26 "y" in base "p", may differ), then the pseudo-sum "z" is represented by O.Z x Z 2 • • - Z L _ X Z L , 

27 where each "Z," is a digit in base "p" given by Z, = ix } + Y t + q) mod p 9 where "mod" 
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represents the modulo function, and C r = 



0 otherwise 



is a carry value from the 



2 "l-lst" digit position, with C, being set to zero. 

3 Using the pseudo-sum function as described above, the generalized Halton sequence that is 

4 used in the system described in the Grabenstein application is generated as follows. If "p" is an 

5 integer, and x 0 is an arbitrary value on the interval from zero to one, inclusive, then the "p fl -adic von 

6 Neumann-Kakutani transformation T p (x) is given by 



1 

T(x) := x® - 
P P P 



(4), 



s 9 



and the generalized Halton sequence x 0 , x„ x 2 is defined recursively as 



(5) 



|§ From equations (4) and (5), it is clear that, for any value for "p," the generalized Halton sequence 

f 1 can provide that a different sequence will be generated for each starting value of "x," that is, for each 

12 x 0 . It will be appreciated that the Halton sequence Hp k as described above is a special case of the 

13 generalized Halton sequence (equations (4) and (5)) for x 0 =0. 

14 The use of a strictly deterministic low-discrepancy sequence can provide a number of 

1 5 advantages over the random or pseudo-random numbers that are used in connection with the Monte 

16 Carlo technique. Unlike the random numbers used in connection with the Monte Carlo technique, 

17 the low discrepancy sequences ensure that the sample points are more evenly distributed over a 

18 respective region or time interval, thereby reducing error in the image which can result from 

19 clumping of such sample points which can occur in the Monte Carlo technique. That can facilitate 

20 the generation of images of improved quality when using the same number of sample points at the 

2 1 same computational cost as in the Monte Carlo technique. 
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1 Summary Of The Invention 

2 The invention provides a new and improved computer graphics system and method for 

3 generating pixel values for pixels in the image using a strictly deterministic methodology for 

4 generating sample points for use in evaluating integrals defining aspects of the image. 

5 In brief summary, the invention provides, in one aspect, a computer graphics system for 

6 generating a pixel value for a pixel in an image, the pixel being representative of a point in a scene 

7 as recorded on an image plane of a simulated camera. The computer graphics system comprises a 
sample point generator and a function evaluator. The sample point generator is configured to 

f*) generate a set of sample points representing at least one simulated element of the simulated camera, 

ffiD the sample points representing elements of, illustratively, for sample points on the image plane, 

U during time interval during which the shutter is open, and on the lens, a Hammersley sequence, and, 

fl for use in global illumination, a scrambled Halton sequence. The function evaluator configured to 

i 3 generate at least one value representing an evaluation of said selected function at one of the sample 

j| points generated by said sample point generator, the value generated by the function evaluator 

tf5 corresponding to the pixel value. 

16 Brief Description Of The Drawings 

17 This invention is pointed out with particularity in the appended claims. The above and 

18 further advantages of this invention may be better understood by referring to the following 

1 9 description taken in conjunction with the accompanying drawings, in which: 

20 FIG. 1 depicts an illustrative computer graphics system constructed in accordance with the 

2 1 inventionn 

22 Detailed Description of an Illustrative Embodiment 
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1 The invention provides an computer graphic system and method for generating pixel values 

2 for pixels in an image of a scene, which makes use of a strictly-deterministic methodology for 

3 generating sample points for use in generating sample values for evaluating the integral or integrals 

4 whose function(s) represent the contributions of the light reflected from the various points in the 

5 scene to the respective pixel value, rather than the random or pseudo-random Monte Carlo 

6 methodology which has been used in the past. The strictly-deterministic methodology ensures a 
1 priori that the sample points will be generally more evenly distributed over the interval or region 

8 over which the integral(s) is (are) to be evaluated in a low-discrepancy manner. 

9 FIG. 1 attached hereto depicts an illustrative computer system 10 that makes use of such a 

n 

P strictly deterministic methodology. With reference to FIG. 1, the computer system 10 in one 

p{ embodiment includes a processor module 1 1 and operator interface elements comprising operator 

H2 input components such as a keyboard 1 2 A and/or a mouse 1 2B (generally identified as operator input 

ffi 

jp element(s) 12) and an operator output element such as a video display device 13. The illustrative 

ft computer system 10 is of the conventional stored-program computer architecture. The processor 

3 

Q module 1 1 includes, for example, one or more processor, memory and mass storage devices, such 

re as disk and/or tape storage elements (not separately shown), which perform processing and storage 

K? operations in connection with digital data provided thereto. The operator input element(s) 12 are 

p provided to permit an operator to input information for processing. The video display device 13 is 

1 9 provided to display output information generated by the processor module 1 1 on a screen 1 4 to the 

20 operator, including data that the operator may input for processing, information that the operator may 

21 input to control processing, as well as information generated during processing. The processor 

22 module 11 generates information for display by the video display device 13 using a so-called 

23 "graphical user interface" ("GUI"), in which information for various applications programs is 

24 displayed using various "windows." Although the computer system 10 is shown as comprising 

25 particular components, such as the keyboard 12A and mouse 12B for receiving input information 

26 from an operator, and a video display device 13 for displaying output information to the operator, 

27 it will be appreciated that the computer system 10 may include a variety of components in addition 

28 to or instead of those depicted in FIG. 1 . 
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1 In addition, the processor module 1 1 includes one or more network ports, generally identified 

2 by reference numeral 1 4, which are connected to communication links which connect the computer 

3 system 10 in a computer network. The network ports enable the computer system 10 to transmit 

4 information to, and receive information from, other computer systems and other devices in the 

5 network. In a typical network organized according to, for example, the client-server paradigm, 

6 certain computer systems in the network are designated as servers, which store data and programs 

7 (generally, "information") for processing by the other, client computer systems, thereby to enable 

8 the client computer systems to conveniently share the information. A client computer system which 

9 needs access to information maintained by a particular server will enable the server to download the 
f§ information to it over the network. After processing the data, the client computer system may also 
33 return the processed data to the server for storage. In addition to computer systems (including the 
MS above-described servers and clients), a network may also include, for example, printers and facsimile 
g§ devices, digital audio or video storage and distribution devices, and the like, which may be shared 
f4 among the various computer systems connected in the network. The communication links 
15 interconnecting the computer systems in the network may, as is conventional, comprise any 
6| convenient information-carrying medium, including wires, optical fibers or other media for carrying 
fj signals among the computer systems. Computer systems transfer information over the network by 
5 means of messages transferred over the communication links, with each message including 
R) information and an identifier identifying the device to receive the message. 

20 It will be helpful to initially provide some background on operations performed by the 

21 computer graphics system in generating an image. Generally, the computer graphic system 

22 generates an image attempts to simulate an image of a scene that would be generated by a camera. 

23 The camera includes a shutter that will be open for a predetermined time T starting at a time t 0 to 

24 allow light from the scene to be directed to an image plane. The camera may also include a lens that 

25 serves to focus light from the scene onto the image plane. The average radiance flux L m n through 

26 a pixel at position (m,n) on an image plane P, which represents the plane of the camera's recording 

27 medium, is determined by 



VTTtL l° +T I L(h(xJ,y)-o)(x,t,y))f mn (x)dydtdx (6) 



A P 
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1 where "A P " refers to the area of the pixel, A L refers to the area of the portion of the lens through 

2 which rays of light pass from the scene to the pixel, and f m n represents a filtering kernel associated 

3 with the pixel. An examination of the integral in equation (6) will reveal that the variables of 

4 integration, "x," "y" and "t," are essentially dummy variables with the variable "y" referring to 

5 integration over the lens area (A L ), the variable "t" referring to integration over time (the time 

6 interval from t 0 to t 0 +T) and the variable "x" referring to integration over the pixel area (A P ). 

7 The value of the integral in equation (6) is approximated by identifying N p sample points x { 

8 in the pixel area, and, for each sample point, shooting N T rays at times t f j in the time interval t 0 to 

9 t 0 +T through the focus into the scene, with each ray spanning N L sample points y ;j k on the lens area 
10 A L . The manner in which subpixel jitter positions x i5 points in time ty and positions on the lens y iJk 

13 are determined will be described below. These three parameters determine the primary ray hitting 

01 

H the scene geometry in h(x i ,t ij ,y ij k ) with the ray direction co^t^y^ k ). In this manner, the value of 
the integral in equation (6) can be approximated as 

fn 

U , =0 i\ T J=0 iy L k=0 

y where "N !t is the total number of rays directed at the pixel. 

1=6 It will be appreciated that rays directed from the scene toward the image can comprise rays 

17 directly from one or more light sources in the scene, as well as rays reflected off surfaces of objects 

18 in the scene. In addition, it will be appreciated that a ray that is reflected off a surface may have 

19 been directed to the surface directly from a light source, or a ray that was reflected off another 

20 surface. For a surface that reflects light rays, a reflection operator T fr is defined that includes a 

21 diffuse portion T fd , a glossy portion T fg and a specular portion T fs , or 

22 T /r = 7>, + 7>, + 7>, (8). 

23 In that case, the Fredholm integral equation L^+TfrL governing light transport can be represented 

24 as 
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2 where transparency has been ignored for the sake of simplicity; transparency is treated in an 

3 analogous manner. The individual terms in equation (9) are 

4 (i) L e represents flux due to a light source; 

5 (ii) Tj_jL e (where Tj _^ - - Tj ) represents direct illumination, that is, flux reflected 

M> off a surface that was provided thereto directly by a light source; the specular component, associated 

ST with the specular portion T fs of the reflection operator, will be treated separately since it is modeled 

ji> as a 8-distribution; 

03 

X§ (iii) T fg ( L - L e ) represents glossy illumination, which is handled by recursive distribution 

%% ray tracing, where, in the recursion, the source illumination has already been accounted for by the 

CP 

fi direct illumination (item (ii) above); 

(iv) T f L represents a specular component, which is handled by recursively using "L" for 

1 3 the reflected ray; 

14 (v) TjTj +jL (where Tj + y =7^ +T f ) represents a caustic component, which is a ray 

1 5 that has been reflected off a glossy or specular surface (reference the T f , f operator) before hitting 

Jg + Js 

16 a diffuse surface (reference the T f operator); this contribution can be approximated by a high 

J d 

1 7 resolution caustic photon map; and 

1 8 (vi) Tj- Tj d L represents ambient light, which is very smooth and is therefore approximated 

1 9 using a low resolution global photon map. 



-11- 



MENT-061 

1 As noted above, the value integral (equation (6)) is approximated by solving equation (7) 

2 making use of sample points x i5 ty and y iJJc , where "x" refers to sample points within area A L of the 

3 respective pixel at location (m,n) in the image plane, "ty" refers to sample points within the time 

4 interval t 0 to t 0 +T during which the shutter is open, and "y y k " refers to sample points on the lens A L . 

5 In accordance with one aspect of the invention, the sample points x { comprise two-dimensional 

6 Hammersley points, which are defined as ^— , O 2 (/) , where 0 < / < N : = (l n ) 2 for a selected 

\N J 

7 integer parameter "n," and <P 2 (i) refers to the radical inverse of "i" in base "two." Generally, the "s" 

8 dimensional Hammersley point set is defined as 



m 



u Hammersley ;{(),...,#- l} -> /' 



(10), 



where I s is the s-dimensional unit cube [0, 1) (that is, an s-dimensional cube that includes "zero," and 



Jj excludes "one"), the number of points "N" in the set is fixed and b„...,b s ., are bases. The bases do 

fl not need to be prime numbers, but they are preferably relatively prime to provide a relatively 

Q 

JL3 uniform distribution. The radical inverse function <D b , in turn, is generally defined as 
14 4>„:iV 0 -> / 



y=0 y=0 



1 5 where I a I is the representation of "i" in integer base "b." The two-dimensional Hammersley 

V J Jj=o 

16 points form a stratified pattern on a 2 n by 2 n grid. Considering the grid as subpixels, the complete 

1 7 subpixel grid underlying the image plane can be filled by simply abutting copies of the grid to each 

18 other. 
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1 Given integer subpixel coordinates (s x ,s y ) the instance "i" and coordinates (x,y) for the sample 

2 point x 4 in the image plane can be determined as follows. Preliminarily, examining 

,y~ i 

3 — J/ , one observes that 



4 (a) each line in the stratified pattern is a shifted copy of another, and 

5 (b) the pattern is symmetric to the line y=x, that is, each column is a shifted copy of another 
iri> column. 

■csx 

s j , s 

1^ Accordingly, given the integer permutation a {k) :~ 2 n $ 2 ( k) for 0<k<2 n , subpixel coordinates 

m 

p (s x ,s y ) are mapped onto strata coordinates (j\k):- is x mod 2 n 9 s y mod 2 n J , an instance number 

0> "i" is computed as 

m 
■™ - 

Sf i = j2" + a(k) (12) 

1 1 and the positions of the jittered subpixel sample points are determined according to 
12 



x, = (s x + d> 2 {k),s y + <t> 2 (j)) 



_ { *(*) iWl (13) 

13 An efficient algorithm for generating the positions of the jittered subpixel sample points x { will be 

14 provided below in connection with Code Segment 1 . A pattern of sample points whose positions 

15 are determined as described above in connection with equations (12) and (13) has an advantage of 

1 6 having much reduced discrepancy over a pattern determined using a Halton sequence or windowed 

1 7 Halton sequence, as described in the aforementioned Grabenstein application, and therefore the 
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approximation described above in connection with equation (7) gives in general a better estimation 
to the value of the integral described above in connection with equation (6). In addition, if "n ?l is 
sufficiently large, sample points in adjacent pixels will have different patterns, reducing the 
likelihood that undesirable artifacts will be generated in the image. 

A ray tree is a collection of paths of light rays that are traced from a point on the simulated 
camera's image plane into the scene. The computer graphics system 10 generates a ray tree by 
recursively following transmission, subsequent reflection and shadow rays using trajectory splitting. 
In accordance with another aspect of the invention, a path is determined by the components of one 
vector of a global generalized scrambled Hammersley point set. Generally, a scrambled 
Hammersley point set reduces or eliminates a problem that can arise in connection with higher- 
dimensioned low-discrepancy sequences since the radical inverse function <D b typically has 

subsequences of b-1 equidistant values spaced by . Although these correlation patterns are 

merely noticeable in the full s-dimensional space, they are undesirable since they are prone to 
aliasing. The computer graphics system 10 attenuates this effect by scrambling, which corresponds 
to application of a permutation to the digits of the b-ary representation used in the radical inversion. 
For a permutation a from a symmetric group S b over integers 0,...,b-l , the scrambled radical inverse 
is defined as 



If the permutation "a" is the identity, the scrambled radical inverse corresponds to the unscrambled 
radical inverse. In one embodiment, the computer graphics system generates the permutation a 
recursively as follows. Starting from the permutation a 2 =(0 5 2) for base b=2, the sequence of 
permutations is defined as follows: 



® b :N 0 xS b -> / 




(14). 
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1 (i) if the base "b" is even, the permutation a b is generated by first taking the values of 2 a h 

2 

2 and appending the values of 2 a h + 1 , and 

2 

3 (ii) if the base "b" is odd, the permutation a b is generated by taking the values of a b . 1? 

b-l 

4 incrementing each value that is greater than or equal to by one, and inserting 

0 the value b-l in the middle. 
%i This recursive procedure results in 

01 o 2 = (0,l) 
3 a 3 = (0,1,2) 
-=? g 4 = (0,2,1,3) 
g o 5 = (0,3,2,1,4) 
Q o 6 = (0,2,4,1,3,5) 
O o 7 = (0,2,5,3,1,4,6) 

f3 o 8 = (0,4,2,6,1,5,3,7).... 

14 The computer graphics system 10 can generate a generalized low-discrepancy point set as 

1 5 follows. It is possible to obtain a low-discrepancy sequence by taking any rational s-dimensional 

16 point "x" as a starting point and determine a successor by applying the corresponding incremental 

1 7 radical inverse function to the components of "x." The result is referred to as the generalized low- 

18 discrepancy point set. This can be applied to both the Halton sequence and the Hammersley 

19 sequence. In the case of the generalized Halton sequence, this can be formalized as 
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1 where the integer vector i 2 , ... , i s ) represents the offsets per component and is fixed in advance 

2 for a generalized sequence. The integer vector can be determined by applying the inverse of the 

3 radical inversion to the starting point "x." A generalized Hammersley sequence can be generated 

4 in an analogous manner. 

5 Returning to trajectory splitting, generally trajectory splitting is the evaluation of a local 

6 integral, which is of small dimension and which makes the actual integrand smoother, which 

7 improves overall convergence. Applying replication, positions of low-discrepancy sample points 

8 are determined that can be used in evaluating the local integral. The low-discrepancy sample points 

9 are shifted by the corresponding elements of the global scrambled Hammersley point set. Since 
O) trajectory splitting can occur multiple times on the same level in the ray tree, branches of the ray tree 
g| are decorrelated in order to avoid artifacts, the decorrelation being accomplished by generalizing the 
f| global scrambled Hammersley point set. 

|| An efficient algorithm for generating a ray tree will be provided below in connection with 

k% Code Segment 2. Generally, in that algorithm, the instance number "i" of the low-discrepancy 

O vector, as determined above in connection with equation (12), and the number "d" of used 

|J) components, which corresponds to the current integral dimension, are added to the data structure that 

# is maintained for the respective ray in the ray tree. The ray tree of a subpixel sample is completely 

pJS specified by the instance number "i." After the dimension has been set to "two," which determines 

1 9 the component of the global Hammersley point set that is to be used next, the primary ray is cast into 

20 the scene to span its ray tree. In determining the deterministic splitting by the components of low 

2 1 discrepancy sample points, the computer graphics system 10 initially allocates the required number 

22 of dimensions "Ad." For example, in simulating glossy scattering, the required number of 

23 dimensions will correspond to "two." Thereafter, the computer graphics system 10 generates 

24 scattering directions from the offset given by the scrambled radical inverses 

25 0 b d (*, * A ■ ■ • , * Wl (i, <? w, ) ' the instances 




-16- 



MENT-061 

where refers to addition modulo "one." Each direction of the "M" replicated rays is determined 
by y { j and enters the next level of the ray tree with d } : = d + A d as the new integral dimension in 

order to use the next elements of the low-discrepancy vector, and i 1 = i + j in order to decorrelate 
subsequent trajectories. Using an infinite sequence of low-discrepancy sample points, the replication 
heuristic is turned into an adaptive consistent sampling arrangement. That is, computer graphics 
system 1 0 can fix the sampling rate AM, compare current and previous estimates every AM samples, 
and, if the estimates differ by less than a predetermined threshold value, terminate sampling. The 
computer graphics system 1 0 can, in turn, determine the threshold value, by importance information, 
that is, how much the local integral contributes to the global integral. 

As noted above, the integral described above in connection with equation (6) is over a finite 
time period T from t 0 to t 0 +T, during which time the shutter of the simulated camera is open. During 
the time period, if an object in the scene moves, the moving object may preferably be depicted in the 
image as blurred, with the extent of blurring being a function of the object's motion and the time 
interval t 0 +T. Generally, motion during the time an image is recorded is linearly approximated by 
motion vectors, in which case the integrand (equation (6)) is relatively smooth over the time the 
shutter is open and is suited for correlated sampling. For a ray instance "i," started at the subpixel 
position x i5 the offset <B 3 (i) into the time interval is generated and the N T -1 subsequent samples 

$ 3 0") + mod 1 for 0<j<N T that is 

1 V <y 



:= t 0 + 



j 



(17). 



Determining sample points in this manner fills the sampling space, resulting in a more rapid 
convergence to the value of the integral (equation (6)). For subsequent trajectory splitting, rays are 
decorrelated by setting the instance i' = i + j. 



In addition to determining the position of the jittered subpixel sample point x ( , and adjusting 
the camera and scene according to the sample point ty for the time, the computer graphics system 
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1 also simulates depth of field. In simulating depth of field, the camera to be simulated is assumed 

2 to be provided with a lens having known optical characteristics and, using geometrical optics, the 

3 subpixel sample point x s is mapped through the lens to yield a mapped point x { \ The lens is sampled 

4 by mapping the dependent samples 



( 



onto the lens area A L using a suitable one of a plurality of known transformations. Thereafter, a ray 
is shot from the sample point on the lens specified by y^ k through the point x- into the scene. The 

j§ offset ( 0 5 (/ + j 9 cr 5 ), 0 7 (/ + j 9 a 1 )) in equation (18) comprise the next components taken from 

Q3 



the generalized scrambled Hammersley point set, which, for trajectory splitting, is displaced by the 



tf> elements 

01 



' k ^ 



of the two-dimensional Hammersley point set. The instance of the ray 



M originating from sample point y tj k is set to "i+j+k" in order to decorrelate further splitting down the 

12 raytree. In equation (18), the scrambled samples 5 (/+ y,<7 5 ),<I> 7 (z" + y,cr 7 )) are used instead 

13 of the unscrambled samples of 5 (i+ j^ 9 0 7 (i+ y)) since in bases "five" and "seven" the up to 

14 five unscrambled samples will lie on a straight line, which will not be the case for the scrambled 

15 samples. 

16 In connection with determination of a value for the direct illumination (7^ _^ L e above), 

1 7 direct illumination is represented as an integral over the surface of the scene dV , which integral is 

1 8 decomposed into a sum of integrals, each over one of the "L" single area light sources in the scene. 

19 The individual integrals in turn are evaluated by dependent sampling, that is 
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{ T fr-fs L e)(y> Z ) = lr L e( X >y) f r {x,y,z)G(x,y)dx 

L 

= il^j)/,^. 2 )^.^ (19) . 

k-\ ^ e ' k 



* I -TT I L e {x j ,y)f r {x j ,y,z)G{x j ,y) 

k=\ 1V1 K 7 = 0 



where suppL e k refers to the surface of the respective "k-th" light source. In evaluating the estimate 
of the integral for the "k-th" light source, for the M k -th query ray, shadow rays determine the fraction 
of visibility of the area light source, since the pont visibility varies much more than the smooth 
shadow effect. For each light source, the emission L e is attenuated by a geometric term G, which 
includes the visibility, and the surface properties are given by a bidirectional distribution function 
f r -f s . These integrals are local integrals in the ray tree yielding the value of one node in the ray tree, 
and can be efficiently evaluated using dependent sampling. In dependent sampling, the query ray 
comes with the next free integral dimension "d" and the instance "i," from which the dependent 
samples are determined in accordance with 



x j 



-fr* 'U'-O® **(>)] < 2o > 



The offset ($> b d [i-> a b d \® b d ^°b d , )) a S a * n * s taken fr° m the corresponding generalized 

scrambled Hammersley point set, which shifts the two-dimensional Hammersley point set 
f ; 



^- 5 0 2 (y)j on the light source. Selecting the sample rate M k = 2" k as a power of two, the 



local minima is obtained for the discrepancy of the Hammersley point set that perfectly stratifies the 
light source. As an alternative, the light source can be sampled using an arbitrarily-chosen number 
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M k -\ 



M k of sample points using 



,0 



as a replication rule. Due to the implicit 



stratification of the positions of the sample points as described above, the local convergence will be 
very smooth. 



The glossy contribution T, ( L - L e J is determined in a manner similar to that described 



above in connection with area light sources (equations (19) and (20), except that a model f g used to 
simulate glossy scattering will be used instead of the bidirectional distribution function f r . In 
determining the glossy contribution, two-dimensional Hammersley points are generated for a fixed 



the current ray tree depth given by the dimension field "d" of the incoming ray. The ray trees 
spanned into the scattering direction are decorrelated by assigning the instance fields i-i+j in a 
manner similar to that done for simulation of motion blur and depth of field, as described above. 
The estimates generated for all rays are averaged by the splitting rate "M" and propagated up the ray 
tree. 

Volumetric effects are typically provided by performing a line integration along respective 
rays from their origins to the nearest surface point in the scene. In providing for a volumetric effect, 

the computer graphics system 10 generates from the ray data a corresponding offset O bd (i) which 

it then uses to shift the M equidistant samples on the unit interval seen as a one-dimensional torus. 
In doing so, the rendering time is reduced in comparison to use of an uncorrelated jittering 
methodology. 

Global illumination includes a class of optical effects, such as indirect illumination, diffuse 
and glossy inter-reflections, caustics and color bleeding, that the computer graphics system 10 
simulates in generating an image of objects in a scene. Simulation of global illumination typically 




splitting rate M and shifted modulo "one" by the offset 



taken from 
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1 involves the evaluation of a rendering equation. For the general form of an illustrative rendering 

2 equation useful in global illumination simulation, namely: 

3 L(x,w) = 4(ic,w)+ w)G{x,x')v(x,?)L{x,w')dA' (21) 

4 it is recognized that the light radiated at a particular point x in a scene is generally the sum of two 

5 components, namely, the amount of light (if any) that is emitted from the point and the amount of 

6 light (if any) that originates from all other points and which is reflected or otherwise scattered from 

J the point x . In equation (21), L(x 9 w) represents the radiance at the point x in the direction 

3 w = ( 0 , (p ) (where ,f 6" represents the angle of direction w relative to a direction orthogonal of the 
surface of the object in the scene containing the point 3c , and "cp" represents the angle of the 

as 

H' _ - _,\ 

M) component of direction w in a plane tangential to the point x ). Similarly, L\x' , w') in the 

3 

|l integral represents the radiance at the point x f in the direction w } - {0\(p*) (where "9'" represents 

fi the angle of direction w } relative to a direction orthogonal of the surface of the object in the scene 

k3 containing the point 3c' , and "cp" 1 represents the angle of the component of direction w } in a plane 

14 tangential to the point x } ), and represents the light, if any, that is emitted from point x } which may 

1 5 be reflected or otherwise scattered from point x . 

1 6 Continuing with equation (2 1 ), L e {x , w) represents the first component of the sum, namely, 

1 7 the radiance due to emission from the point x in the direction w , and the integral over the sphere 

18 S' represents the second component, namely, the radiance due to scattering of light at point x . 

19 f{x , w y -> w) is a bidirectional scattering distribution function which describes how much of the 

20 light coming from direction w' is reflected, refracted or otherwise scattered in the direction w , and 

-21- 




I 

MENT-061 

1 is generally the sum of a diffuse component, a glossy component and a specular component. In 

2 equation (21), the function G(x, x') is a geometric term 



, . COS0 cos#' .... 

G(x,x>) = -pr^jr- (22) 

be - x 



4 where 9 and 9'are angles relative to the normals of the respective surfaces at points x and x' , 

5 respectively. Further in equation (21), V(x, x') is a visibility function which equals the value one 
^ if the point x' is visible from the point x and zero if the point x' is not visible from the point x . 

02 The computer graphics system 1 0 makes use of global illumination specifically in connection 

S with determination of the diffuse component T fd T fj L and the caustic component T fd T fgJrf L 



pl using a photon map technique. Generally, a photon map is constructed by simulating the emission 

10 of photons by light source(s) in the scene and tracing the path of each of the photons. For each 

y simulated photon that strikes a surface of an object in the scene, information concerning the 

tj simulated photon is stored in a data structure referred to as a photon map, including, for example, 

If the simulated photon's color, position and direction angle. Thereafter a Russian roulette operation 

k% is performed to determine the photon's state, that is, whether the simulated photon will be deemed 

1 5 to have been absorbed or reflected by the surface. If a simulated photon is deemed to have been 

16 reflected by the surface, the simulated photon's direction is determined using, for example, a 

1 7 bidirectional reflectance distribution function ("BRDF"). If the reflected simulated photon strikes 

18 another surface, these operations will be repeated (reference the aforementioned Grabenstein 

19 application). The data structure in which information for the various simulated photons is stored 

20 may have any convenient form; typically k-dimensional trees, for "k" an integer" are used. After the 

2 1 photon map has been generated, it can be used in rendering the respective components of the image. 

22 In generating a photon map, the computer graphics system 1 0 simulates photons traj ectories, 

23 thus avoiding the necessity of discretizing the kernel of the underlying integral equation. The 

24 interactions of the photons with the scene, as described above, are stored and used for density 
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1 estimation. The computer graphics system 1 0 makes use of a scrambled Halton sequence, which has 

2 better discrepancy properties in higher dimensions than does an unscrambled Halton sequence. The 

3 scrambled Halton sequence also has the benefit, over a random sequence, that the approximation 

4 error decreases more smoothly, which will allow for use of an adaptive termination scheme during 

5 generation of the estimate of the integral. In addition, since the scrambled Halton sequence is strictly 

6 deterministic, generation of estimates can readily be parallelized by assigning certain segments of 

7 the low-discrepancy sequence to ones of a plurality of processors, which can operate on portions of 

8 the computation independently and in parallel. Since usually the space in which photons will be shot 

9 by selecting directions will be much larger than the area of the light sources from which the photons 
10 were initially shot, it is advantageous to make use of components of smaller discrepancy, for 
ljj example, G> 2 or <D 3 (where, as above, O b refers to the radical inverse function for base "b"), for 
|2 angular scattering and components of higher discrepancy, for example, <I> 5 or G> 7 for area sampling, 
iB which will facilitate filling the space more uniformly. 

fn 

ijf The computer graphics system 1 0 estimates the radiance from the photons in accordance with 

Q 1 L,(X 9 0) * — Ysfr{®i> X >®)®i (23) 

A ieB k (x) 

Mi where, in equation (23), 0 { represents the energy of the respective "i-th" photon, C0j is the direction 

1 7 of incidence of the "i-th photon, "B k (x) represents the set of the "k" nearest photons around the point 

18 "x," and "A" represents an area around point "x" that includes the photons in the set B k (x). The 

1 9 computer graphics system 1 0 makes use of an unbiased but consistent estimator is used for * * , which 

20 is determined as follows. Taking the radius r(B k (x)) of the query ball (* * what is this?* *) a tangential 

21 disk D of radius r(B k (x)) centered on the point x is divided into M equal-sized subdomains D i9 that 

22 is 



23 



u^" 1 D; = D and D,n D, * 0 for it j, where D\= " ' \ (24) ' 



|Z3| nr 2 (B h (x)) 
M M 



24 The set 
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1 P = {A|^n{xJ D |/e^(x)}*oJ (25) 

2 contains all the subdomains D { that contain a point x^ on the disk, which is the position of the "i-th" 

3 photon proj ected onto the plane defined by the disk D along its direction of incidence co- . Preferably, 

4 the number M of subdomains will be on the order of *Jk and the angular subdivision will be finer 

5 than the radial subdivision in order to capture geometry borders. The actual area A is then 

6 determined by 

7 A = *r 2 (**M)5 (26) ' 

M| Determining the actual coverage of the disk D by photons significantly improves the radiance 

4- 

Of estimate (equation (23)) in corners and on borders, where the area obtained by the standard estimate 

01 

hB n r 2 {B k {x)) would be too small, which would be the case at corners, or too large, which would 

y be the case at borders. In order to avoid blurring of sharp contours of caustics and shadows, the 

tf computer graphics system 10 sets the radiance estimate L to black if all domains D { that touch x do 

jj not contain any photons. 

14 It will be appreciated that, in regions of high photon density, the "k" nearest photons may 

15 lead to a radius r(B k (x) that is nearly zero, which can cause an over-modulation of the estimate. 

16 Over-modulation can be avoided by selecting a minimum radius r min , which will be used if r(B k (x)) 

17 is less than r min . In that case, instead of equation (23), the estimate is generated in accordance with 

L r (x>a) = T L^ifX^i^r) (27) 

19 assuming each photon is started with J/^ of the total flux The estimator in equation (27) 

20 provides an estimate for the mean flux of the "k" photons if r(B k (x))<r min . 
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1 The global photon map is generally rather coarse and, as a result, subpixel samples can result 

2 in identical photon map queries. As a result, the direct visualization of the global photon map is 

3 blurry and it is advantageous to perform a smoothing operation in connection therewith In 

4 performing such an operation, the computer graphics system 1 0 performs a local pass integration that 

5 removes artifacts of the direct visualization. Accordingly, the computer graphics system 10 

6 generates an approximation for the diffuse illumination term 7J. 7J. L as 

7 T fJf< L * { T f£fc) = \ s , (x) fA^)L r {h{x,m))cosedco 

f (x) M ~ l — l ( ( \\\ ^ 28 ^ 

O * ~ A ^f Z L r y\ x > 0 (a««in ju~ x 2nu i2 Jjj 

W 

% with the integral over the hemisphere S (x) of incident directions aligned by the surface normal in 

ffl "x" being evaluated using importance sampling. The computer graphics system 10 stratifies the 

5i sample points on a two-dimensional grid by applying dependent trajectory splitting with the 

£L Hammersley sequence and thereafter applies irradiance interpolation. Instead of storing the incident 

W flux Oj of the respective photons, the computer graphics system 10 stores their reflected diffuse 

l| power fiCx^Oj with the respective photons in the photon map, which allows for a more exact 

P approximation than can be obtained by only sampling the diffuse BRDF in the hit points of the final 

15 gather rays. In addition, the BRDF evaluation is needed only once per photon, saving the 

16 evaluations during the final gathering. Instead of sampling the full grid, the computer graphics 

1 7 system 1 0 uses adaptive sampling, in which refinement is triggered by contrast, distance traveled by 

18 the final gather rays in order to more evenly sample the projected solid angle, and the number of 

19 photons that are incident form the portion of the projected hemisphere. The computer graphics 

20 system fills in positions in the grid that are not sampled by interpolation. The resulting image matrix 

2 1 of the projected hemisphere is median filtered in order to remove weak singularities, after which the 

22 approximation is generated . The computer graphics system 10 performs the same operation in 

23 connection with hemispherical sky illumination. 

24 The computer graphics system 10 processes final gather rays that strike objects that do not 

25 cause caustics, such as plane glass windows, by recursive ray tracing. If the hit point of a final 
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gather ray is closer to its origin than a predetermined threshold, the computer graphics system 10 
also performs recursive ray tracing. This reduces the likelihood that blurry artifacts will appear in 
corners, which might otherwise occur since for close hit points the same photons would be collected, 
which can indirectly expose the blurry structure of the global photon map. 

Generally, photon maps have been taken as a snapshot at one point in time, and thus were 
unsuitable in connection with rendering of motion blur. Following the observation that averaging 
the result of a plurality of photon maps is generally similar to querying one photon map with the 
total number of photons from all of the plurality of photon maps, the computer graphics system 10 
generates N T photon maps, where N T is determined as described above, at points in time 

* = t + (29) ' 

for 0<b<N T . During rendering, the computer graphics system 10 uses the photon map with the 



smallest time difference 



in connection with rendering for the time sample point ty. 



The invention provides a number of advantages. In particular, the invention provides a 
computer graphics system that makes use of strictly deterministic distributed ray tracing based on 
low-discrepancy sampling and dependent trajectory splitting in connection with rendering of an 
image of a scene. Generally, strictly deterministic distributed ray tracing based on low-discrepancy 
sampling and dependent trajectory splitting is simpler to implement than an implementation based 
on random or pseudo-random numbers. Due to the properties of the radical inversion function, 
stratification of sample points is intrinsic and does not need to be considered independently of the 
generation of the positions of the sample points. In addition, since the methodology is strictly 
deterministic, it can be readily parallelized by dividing the image into a plurality of tasks, which can 
be executed by a plurality of processors in parallel. There is no need to take a step of ensuring that 
positions of sample points are not correlated, which is generally necessary if a methodology based 
on random or pseudo-random numbers is to be implemented for processing in parallel. 
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1 Generally, a computer graphics system that makes use of low-discrepancy sampling in 

2 determination of sample points will perform better than a computer graphics system that makes use 

3 of random or pseudo-random sampling, but the performance may degrade to that of a system that 

4 makes use of random or pseudo-random sampling in higher dimensions. By providing that the 

5 computer graphics system performs dependent splitting by replication, the superior convergence of 

6 low-dimensional low-discrepancy sampling can be exploited with the effect that the overall 

7 integrand becomes smoother, resulting in better convergence than with stratified random or pseudo- 

8 random sampling. Since the computer graphics system also makes use of dependent trajectory 

9 sampling by means of infinite low discrepancy sequences, consistent adaptive sampling of, for 
10 example, light sources, can also be performed. 

IS In addition, it will be appreciated that, although the computer graphics system has been 

fif described as making use of sample points generated using generalized scrambled and/or unscrambled 

+= 

W Hammersley and Halton sequences, it will be appreciated that generally any (t,m,s)-net or (t,s)- 

12 sequence can be used. 

Eg Code Segment 1 

CP 

ft The following is a code fragment in the C++ programming language for generating the 

O positions of the jittered subpixel sample points x { 

18 unsigned short Period, *Sigma; 

1 9 void InitSigma(int n) 

20 { 

21 unsigned short Inverse, Digit, Bits; 

22 Period = 1 « n; 

23 Sigma = new unsigned short [Period] ; 

24 for (unsigned short i = 0; i < Period; i++) 

25 { 

26 Digit = Period 

27 Inverse = 0; 

28 for (bits = i; bits; bits »= 1) 
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1 { 

2 Digit »=1; 

3 if (Bits & 1) 

4 inverse += Digit; 

5 } 

6 Sigma[i] = Inverse; 

7 } 

8 } 

9 void SampleSubpixel(unsigned int *i, double *x, double *y, int s x , int s y ) 
10 { 

1? intj = s x & (Period -1); 

IflS int k = s y & (Period - 1 ); 

fjj *i = j * Period + Sigma[k] 

ljj *x = (double) s x + (double) Sigma[k] / (double) Period; 

i *y = (double) s y + (double) Sigmajj] / (double) Period; 

fe } 

fn 
■=j • 

f£ Code Segment 2 

y The following is a.code fragment in the C++ programming language for generating a ray tree 

1 9 class Ray 

20 { 

2 1 int i; //current instance of low discrepancy vector 

22 int d; //current integral dimension in ray tree 
23 

24 } 

25 void Shade (Ray& Ray) 

26 { 

27 Ray next_ray; 

28 int i = ray.i; 

29 int d = ray.d 
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1 

2 for(intj = 0;j <M;j++) 

3 { 
4 

5 // ray set up for recursion 



6 




7 next ray = SetUpRay(y); // determine ray parameters by y 

8 next ray.i = i + j; // decorrelation by generalization 
P nextray.d = d + Ad; //dimension allocation 

l=§j Shade(next_ray); 

lffi 

m } 

f } 

H It will be appreciated that a system in accordance with the invention can be constructed in 

W\ whole or in part from special purpose hardware or a general purpose computer system, or any 

fpi combination thereof, any portion of which may be controlled by a suitable program. Any program 

Q may in whole or in part comprise part of or be stored on the system in a conventional manner, or it 

18 may in whole or in part be provided in to the system over a network or other mechanism for 

1 9 transferring information in a conventional manner. In addition, it will be appreciated that the system 

20 may be operated and/or otherwise controlled by means of information provided by an operator using 

21 operator input elements (not shown) which may be connected directly to the system or which may 

22 transfer the information to the system over a network or other mechanism for transferring 

23 information in a conventional manner. 

24 The foregoing description has been limited to a specific embodiment of this invention. It will 

25 be apparent, however, that various variations and modifications may be made to the invention, with 

26 the attainment of some or all of the advantages of the invention. It is the object of the appended 

27 claims to cover these and such other variations and modifications as come within the true spirit and 

28 scope of the invention. 
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What is claimed as new and desired to be secured by Letters Patent of the United States is: 
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